GATK之HaplotypeCaller

GATK的主要功能其实就是识别变异位点，其他功能都是锦上添花。所以这一次学习GATK寻找变异位点的工具。
在GATK的文档中，与变异位点识别相关的有9个工具，分别是：

Name	Summary
	Apply a score cutoff to filter variants based on a recalibration table
	Calculate genotype posterior likelihoods given panel data
	Simple Bayesian genotyper used in the original GATK paper
	Perform joint genotyping on gVCF files produced by HaplotypeCaller
	Call germline SNPs and indels via local re-assembly of haplotypes
	Call somatic SNPs and indels via local re-assembly of haplotypes
	Regenotypes the variants from a VCF containing PLs or GLs.
	Call SNPs and indels on a per-locus basis
	Build a recalibration model to score variant quality for filtering purposes

稍微看了一下GATK最佳实践文档，发现目前用的最多的是HaplotypeCaller，准确率高，但是比较吃配置，所以运行时间会比较久。不过由于HaplotypeCaller的工作原理，直接省去了BQSR和indel realignment步骤，所以对于一个variant calling流程而言，可以直接比对，去重复后运行HaplotypeCaller。

简介

功能： 通过局部单倍型重组装寻找体细胞的SNP和indel
分类： 变异位点探索功能

HaplotypeCaller，简称HC，能过通过对活跃区域（也就是与参考基因组不同处较多的区域）局部重组装，同时寻找SNP和INDEL。换句话说，当HC看到一个地方好活跃呀，他就不管之前的比对结果了，直接对这个地方进行重新组装，所以就比传统的基于位置（position-based)的工具，如前任UnifiedGenotyper好用多了。

HC可以处理非二倍体物种和混池实验数据，但是不推荐用于癌症或者体细胞。
HC还可以正确处理可变剪切，所以也能用于RNA-Seq数据。

HC有一个GVCF模式，产生各个样本的中间基因组文件(gVCF),然后用于多样本的联合基因定型(joint genotyping)，效率很高呢。

工作原理

HaplotypeCaller的核心操作就是四步：

1. 寻找活跃区域，就是和参考基因组不同部分较多的区域

2. 通过对该区域进行局部重组装，确定单倍型（haplotypes）。就是这一步可以省去indel realignment

3. 在给定的read数据下，计算单倍型的可能性。

4. 分配样本的基因型

使用说明

输入：BAM文件
输出：原始的VCF文件或者gVCF文件，未过滤SNP和INDEL。一般而言，在进行下游分析，需要对这些数据要么手动要么VQSR过滤。

案例:
DNA-seq variant-calling

  java -jar GenomeAnalysisTK.jar \     -R reference.fasta \     -T HaplotypeCaller \     -I sample1.bam [-I sample2.bam ...] \     [--dbsnp dbSNP.vcf] \     [-stand_call_conf 30] \     [-L targets.interval_list] \     -o output.raw.snps.indels.vcf

RNA-seq variant-calling

  java -jar GenomeAnalysisTK.jar \     -R reference.fasta \     -T HaplotypeCaller \     -I sample1.bam \     [--dbsnp dbSNP.vcf] \     -stand_call_conf 20 \     -o output.raw.snps.indels.vcf

其他感觉比较使用的参数：

参数名	默认值	概要
—genotyping_mode/-gt_mode	DISCOVERY	如何处理识别的等位基因
—min_base_quality_score/ -mbq	10	最低碱基质量
—sample_ploidy/-ploidy	2	样本是几倍体？
—standard_min_confidence_threshold_for_calling/-stand_call_conf	10.0	确定variant的最低phred-scaled 置信阈值
—annotation/-A	-	variant注释内容
—pcr_indel_model/-pcrModel	CONSERVATIVE	PCR indel模式

注：对于我做mapping-by-sequencing而言，需要结果有ref和alt碱基的支持数，所以选项-A一定要跟上StrandAlleleCountsBySample。

这是我用于MBS的流程的GATK部分

############################## # Vaiant Discovery with GATK # ##############################mkdir -p ${wd_dir}/analysis/variant_gatk recal_files=$(ls *recal_reads.bam) for file in $recal_files do    java -Xmx16g -jar $gatk -T HaplotypeCaller \    -R $reference  -I $file --genotyping_mode DISCOVERY \    --standard_min_confidence_threshold_for_calling 10 \    -A StrandAlleleCountsBySample \    -o ../variant_gatk/${file%%.*}_raw_variants.vcf done

更多参数详细说明见，可以阅读原文。

写在最后

GATK best practice写在2015年，现在是2017年，GATK的版本到3.7，下一个版本就是4.0了，所以很多工具都不能想当然的用了。比如说之前写的BQSR，有一篇文献就说已经没有必要了，文献名见下：
《 impact of post-alignment processing in variant discovery from whole exome data>

而这个也是我发完BQSR，留言的大神告诉我的。后来我去biostar搜索，也发现有人说BQSR和VQSR都不是必须的（我正在在问回答者原文出处）

技术的发展越来越快，对于一个刚入门的bioinformatician，很多东西都要从头开始学习，压力非常大。所以我非常感谢一个论坛—技能树，非常感谢一本书—biostar handbook.还有我们生信媛的小伙伴们。

没事多去逛biostar，seqanswers和生信技能树，开拓自己的眼界。也希望那些老司机们不要吝啬自己的知识，能够带带我们。

最后，让我们“一起学习，共同进步”